Making a map with points in R with leaflet

We’re going to be doing a lot of the stuff from the RStudio leaflet tutorial. This walkthrough is a variation of Amelia McNamara’s version (you’ll meet her later this summer).

Load leaflet first:

library(leaflet)

Let’s start with a simple map

m <- leaflet() %>%
  setView(-72.518978, 42.381050, zoom = 6) %>%
  addTiles() %>%  # Load the default OpenStreetMap tile set (images that make up the map)
  addMarkers(lng=-72.518978, lat=42.381050, popup="DSDP!")

m

Note that the syntax here is similar to that of ggplot2. Please consult the Leaflet documentation for more details.

Okay, let’s do something with actual data

Lets look at storm data from the NOAA. It comes in a few files that we need to join together in order to use. For convenience, I’ve put them on the workshop website. You can download them manually, or you can use the getURL() function from RCurl to grab them from the web as I do below (they’re big, so it takes a sec):

library(readr)
library(dplyr)
library(RCurl)
## Loading required package: bitops
## 
## Attaching package: 'RCurl'
## The following object is masked from 'package:tidyr':
## 
##     complete
# One file contains information about the storm
stormdetails <- getURL("https://jcrouser.github.io/MassMutual-DataVis/datasets/StormEvents_details-ftp_v1.0_d2016_c20160810.csv")
stormdetails <- read_csv(stormdetails)

# The other contains actual location data
stormlocs <- getURL("https://jcrouser.github.io/MassMutual-DataVis/datasets/StormEvents_locations-ftp_v1.0_d2016_c20160810.csv")
stormlocs <- read_csv(stormlocs)

# We'll want to use them together, so we'll use a join
storms <- stormlocs %>%
  left_join(stormdetails, by="EVENT_ID")

Don’t worry if you don’t really understand the syntax above. Next week, Ben Baumer will teach you all about SQL and joins.

Mapping lightning strikes

Let’s pull out the lightning strikes and map them as part of a dplyr chain

lightning_map <- storms %>%               # Start with the storm data
  filter(EVENT_TYPE == "Lightning") %>%   # Filter down to just the lightning events
  leaflet() %>%                           # Pipe into a leaflet map
    addMarkers(~LONGITUDE, ~LATITUDE) %>%   # Add markers at each LONGITUDE/LATITUDE pair
    addProviderTiles("Stamen.Toner")        # We'll use black/white tiles for dramatic effect

lightning_map

Challenge:

  • Find another storm type to map
  • Bonus– add popups!

One approach

mtw_map <- storms %>%
  filter(EVENT_TYPE == "Marine Thunderstorm Wind") %>%
  leaflet() %>%
    addProviderTiles("Stamen.Toner") %>% 
    addMarkers(~LONGITUDE, ~LATITUDE, popup = ~EVENT_NARRATIVE)

mtw_map

Polygons

First let’s grab some data (and do a little conversion)

tornados <- storms %>%
  filter(EVENT_TYPE=="Tornado") %>%
  mutate(DAMAGE_PROPERTY = as.numeric(sub("K", "", DAMAGE_PROPERTY, fixed = TRUE)))
## Warning in eval(substitute(expr), envir, enclos): NAs introduced by
## coercion

Then let’s start with something easy(ish) – circles

m <- leaflet(data = tornados) %>%
  addProviderTiles("Stamen.Toner") %>% 
  addCircles(~LONGITUDE, ~LATITUDE, 
             weight = 1, 
             radius = ~DAMAGE_PROPERTY*100, # Map the radius of the circle to amount of damage
             popup = ~EVENT_NARRATIVE)      # Include details about the tornado
m

Polgyons come in shapefiles

Most boundaries (state, national, etc) are provided in terms of polygons. Major mapping software ArcGIS, from ESRI, has essentially set the standard formats. There are many files with different extensions: .prj (the projection), .shp (the shapefile), .cpg (??), .dbf (??), .shx (??).

You need special software or packages to work with shapefiles.

State shapefiles

I got these from the Census. You can choose the resolution.

If you want, the zipfile of the shapes I used is here.

We’re going to use the maptools package to deal with shapefiles.

library(maptools)
## Loading required package: sp
## Checking rgeos availability: TRUE
# Remember when we talked about map projections?
crswgs84 = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

# Now we can load the shapefile using the correct projection
states = readShapePoly("datasets/cb_2015_us_state_500k/cb_2015_us_state_500k.shp",
                       proj4string = crswgs84,
                       verbose = TRUE)
## Warning: use rgdal::readOGR or sf::st_read
## Shapefile type: Polygon, (5), # of Shapes: 56

In RStudio, you can click on the states object to see what a Large SpatialPolygonsDataFrame looks like. It contains both @data and @polygons.

Exploring shapefiles

Let’s start by mapping some boring, internal data from shapefile about the amount of water per state

states %>%
  leaflet() %>%
  setView(-95.976807, 40.829587, zoom = 3) %>%
  addProviderTiles("Stamen.Toner") %>%
  addPolygons(stroke = FALSE, 
              fillOpacity = 0.5, 
              smoothFactor = 0.5, 
              color = ~colorQuantile("BrBG", states$AWATER)(AWATER)
  )

Challenge

Put it all together! See if you can make a map of your own data, and embed it within a dashboard. Bonus: can you make them interact with one another?